Close packing density of polydisperse hard spheres 
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The most efficient way to pack equally sized spheres isotropically in 3D is known as the random 
close packed state, which provides a starting point for many approximations in physics and engi- 
neering. However, the particle size distribution of a real granular material is never monodisperse. 
Here we present a simple but accurate approximation for the random close packing density of hard 
spheres of any size distribution, based upon a mapping onto a one-dimensional problem. To test this 
theory we performed extensive simulations for mixtures of elastic spheres with hydrodynamic fric- 
tion. The simulations show a general (but weak) dependence of the final (essentially hard sphere) 
packing density on fluid viscosity and on particle size, but this can be eliminated by choosing a 
specific relation between mass and particle size, making the random close packed volume fraction 
well-defined. Our theory agrees well with the simulations for bidisperse, tridisperse and log-normal 
distributions, and correctly reproduces the exact limits for large size ratios. 



I. INTRODUCTION 

Granular materials such as sediments and powders 
are widespread in nature and industrial contexts, and 
treating the grains as hard spheres is often a useful 
first approximation. In these systems, the manner in 
which the grains pack together has profound influence 
on properties. Of these, random close packing[l| is most 
likely to be encountered in tapped and consolidated sys- 
tems, although other possibilities, such as random loose 
packing 2], and various crystalline arrangements (whose 
existence is very sensitive to the form of the size distri- 
bution and method of creation 0, are also possible. 

Even though the precise nature (and for monodisperse 
spheres, even the existence [||]) of the random close packed 
state remains the subject of ongoing research [||, it pro- 
vides a starting point for many approximations in both 
physics 0,Q and engineering, and has great practical im- 
portance not only for the prediction of the density of 
granular materials, but also other properties. For exam- 
ple, the viscosity of dense dispersions will diverge at this 
poin t [9| and it is related to the permeability in packed 
beds[10l|. 

One of the insights that have come forward from 
simulations [1] is that the dense random packing density 
of hard spheres depends upon the (shear) friction coef- 
ficient, if the particles only lose energy by inelastic col- 
lisions. Further, the jamming density of a hard sphere 
system depends upon the initial state, and on the par- 
ticular pathway chosen to cool down the system. In gen- 
eral, however, dissipative interactions play a role not just 
at contact. Granular particles suspended in a viscous 
medium also dissipate energy via long-range hydrody- 
namic interactions. Hence we anticipate that the dense 
random packing also depends upon solvent viscosity and 
on the range of the (hydrodynamic) friction. This is the 
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FIG. 1: Close packed configuration of spheres from a log- 
normal distribution. The spread in the logarithm of radius 
is a = 0.6. Only spheres with centres lying in one periodic 
image of the simulation cell are shown. 



first problem that we wish to address. To this end we 
developed a new simulation method that includes these 
effects. Using this method we not only find that the 
dense random packing depends on fluid viscosity, but - 
quite unexpectedly - also on particle size and mass. By 
analysing the various time scales in the problem we ob- 
tain a way to eliminate this dependence, which sheds new 
light upon the nature of the dense random packed state. 

In practical cases the particle size distribution of a 
real granular material, or mixture of materials, is never 
monodisperse. Also for such polydisperse problems mod- 
elling techniques [111, [l2j have been used to calculate max- 
imum packing fractions of spheres with a distribution of 
sizes. A typical example of a polydisperse system in a 
close packed state is shown in Fig. [1] 

The second problem that we wish to address is that 
these simulation methods are quite time consuming, and 
therefore their applicability is limited. It would be desir- 
able to have an analytic expression for the close packing 
density, or a fast approximation, but progress in this di- 
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rection has not been rapid. Ouchiyama and Tanaka[l3j 
have presented a theory based upon the volume occupied 
by a sphere in contact with other spheres of the mean 
diameter, but their results are at best qualitative, and 
the reasoning behind the method is not simple enough to 
suggest obvious improvements. Song et a/.[6( presented a 
theory for the packing of monodisperse spheres, but the 
generalisation to an arbitrary size distribution is not ob- 
vious. Recently Biazzo et a/.jlj] presented a theory for 
binary mixtures, but like the Ouchiyama- Tanaka theory, 
it violates the exact upper limit for large fractions of big 
spheres that is given in Eq. ([3]) below. Thus, a com- 
prehensive theory to predict the random close packing 
density of an arbitrary sphere mixture is still lacking. 

We formulate such a theory, which is presented here 
in Section II. Next, we define our simulation method in 
Section III, and we present the combined theoretical and 
simulation results in Section IV. Conclusions are formu- 
lated in Section V. 



II. THEORY 

Here we propose an approximate solution to the prob- 
lem of polydisperse packing density, obtained by ab- 
stracting what we believe to be essential features of the 
physics and geometry of packing. The fundamental prob- 
lem we wish to solve is as follows: suppose we have a 
normalised distribution P 3 d{D) of sphere diameters de- 
fined so that P-$d{D)&D equals the number fraction of 
spheres with diameters in the range (D, D + dD) present 
in the system. Then we ask what is the functional 
F '■ Psd(D) l— * '/'max that maps the size distribution onto 
the random close packed volume fraction? 

In order to construct this functional, we begin by 
mapping the 3D sphere packing problem onto a pack- 
ing problem of rods in ID. The corresponding ID dis- 
tribution P\£){L)AL gives the number fraction of rods 
with rod length in the range [L, L + dL) present in the 
system. To do this we imagine a large random, but non- 
overlapping, arrangement of spheres in 3D with size dis- 
tribution P 3 d(D). This need not be close packed for the 
argument that follows. Now imagine drawing a straight 
line through this distribution, and counting each portion 
of the line which lies within a sphere as a rod (see Fig. 
[3]) . The resulting distribution of rod lengths is then given 
by 



P 1D (L) = 2L 



/„" P 3D (D)D^dD' 



(1) 



If rods of length L, are arranged on a line of length A 
(with periodic boundary conditions) , then the rod length 
fraction is clearly given by rp = A J^Z,-. This equals 
the volume fraction 4> if there is a corresponding system 
of spheres in 3D. 

Let us assume that it is possible to map the closest ran- 
dom packing of spheres in 3D onto a problem of packing 




FIG. 2: How to map a 3D sphere distribution onto a rod 
distribution. A straight line through a random arrangement 
of spheres defines a set of rods. The probability that a sphere 
is hit by the line is proportional to its cross sectional area 
and so proportional to D 2 Pzd(D). The probability that such 
a hit produces a rod of length L is 2L/D 2 (for L < D), and 
Eq. (fl| in the text follows. 



the above collection of rods on a line, where we search 
over all orderings of rods as well as their separations. 
This mapping will be achieved through an effective po- 
tential between the rods, which must have the following 
properties: 

1. It should lead to a maximum packing fraction ipmax 
which is unchanged if all the rods (or spheres) are mag- 
nified by an equal amount; 

2. The potential should be 'hard', in that it is either 
zero or infinite; 

3. The interaction potential between large rods should 
reach through small rods. This will allow very small rods 
to 'rattle around' in the gaps between the large rods, so 
that at high weight fractions of large rods the latter can 
form a stress-bearing network; 

4. The interaction range for mixtures of very unequal 
rods should be determined by the size of the smallest rod, 
so that small rods can form a dense randomly packed 
system in between the large rods, without leaving large 
gaps. 

The true interaction potential between the rods will 
be both many-body and highly complicated, capturing 
topological aspects of 3D space. However, we suggest 
that the following pair potential is the simplest expres- 
sion which honours the four requirements listed: 

If two rods i and j have a gap h between their nearest 
approaching ends, then we introduce the potential 



V(h) = 



oo 




if h < min(fLi,fLj) 
if h > min(fLi,fLj) 



(2) 



where / > is a 'free volume' parameter. If this potential 
applies between any pair of rods, regardless of interven- 
ing rods in the gap between them, then it will satisfy 
requirement 3. The other requirements follow naturally. 
Rather than taking the minimum of the two lengths one 
could introduce a more complicated function of and 
Lj, but the present choice appears to be the simplest to 
capture the physics. In the remainder of this article / 
will be used as a fit parameter; it is the only parameter 
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in our theory, and can be chosen by requiring that the 
theory reproduces the random close packing of monodis- 
perse systems. 

For each ordering of the rods on the line, there is a 
shortest line which can accommodate them without in- 
curring an infinite potential energy, and this leads to a 
close packing fraction for this ordering, which is simply 
ip = A -1 ^2 Li- The maximum packing fraction ip max is 
the maximum value attained by tp over all possible order- 
ings of the rods. If all the L's are equal, V'max = (l+.f)~\ 
and any rod polydispersity (which will always be present 
if we use Eq. ((T|)) will increase V'max- 

If we imagine inserting rods one at a time to form a 
packing, while increasing A if necessary, then Eq. ([2|) 
constitutes a two-body potential between the inserted 
rod and all the rods currently in the packing. However, 
if rods are inserted in decreasing order of size, then the 
special choice of this potential means it only depends 
upon the newly inserted rod, and the packing further 
away is not disrupted by this process. To insert the new 
rod with a minimum increase of line length A we need to 
identify the biggest gap. Therefore the following 'greedy 
algorithm' [15| may be used to find ?/> max for arbitrary 
{Li}: 

(a) The set of lengths {Li} is labelled such that L\ > 
L2 > . • • > Ln ■ These will be inserted in decreasing or- 
der of lengths into the growing optimal packing, starting 
with L\. 

(b) Throughout the algorithm, we maintain a set of 
gaps {gi}, equal in number to the number of rods we have 
inserted into the packing. At the start, when we have 
only one rod, this set contains one element g\ = fL\. 

(c) In order to insert rod j, we identify g max , the largest 
gap in the set of gaps, and we remove it from the set. 
We then add two new gaps to the set, namely fLj and 
max[g max - (1 + f)Lj, fLj]. 

This process implicitly increases the line length A, if 
that is needed to accommodate the new rod. 

Our final approximation consists of choosing a large 
number of rods from the distribution Piu(L), and pack- 
ing them by the greedy ID packing algorithm to obtain 
our estimate ijj max for </> max , the close packed volume frac- 
tion of the sphere distribution P$d- Since this is essen- 
tially a sorting routine, the predictions of this algorithm 
take much less time (ca 0.3 seconds on a modern desktop 
computer) than explicit 3D simulation (1 to 30 hours). 
When the rod lengths are chosen at equidistant values of 
the cumulative ID distribution, 2000 rods are sufficient 
for 5 decimal places accuracy. Thus, for N rods we choose 
rod length L % such that f^ 3 P 1D {l)dl = (N - i + 1/2) /N. 
We use N = 20000 rods. 

One useful property of this procedure is that it cor- 
rectly reproduces the exact solution for bidisperse spheres 
with infinite size ratio. This limit is given by 

, • ( 4>RCP 4>RCp\ ,„k 

V>max = mm r, , 3 

\1 - U>(1 - (j>RCp) w ) 

where 4>rcp is the maximum packing fraction for a 



monodisperse system, and w is the mass fraction of 
large spheres on the total particle volume, so w — 
</'iarge/(0iargc + Ismail)- Numerical solution of the the- 
ory for size ratios down to 1:1000 shows minor (~ 1%) 
deviations from the exact limit, for w values very close 
to the cusp, w = 1/(2 — 4>rcp)- 

The description above is complete, except that we 
need to specify the parameter /, which should be chosen 
such that the predicted maximum packing for monodis- 
perse spheres is the correct random close packing value 
4>RCP- For monodisperse spheres we have P\d(L) = 
2LDq 2 6(Dq — L) (where 9 is the Heaviside step func- 
tion), and we find that a value of / = 0.7654 leads to 
a packing fraction of approximately 0.6435. This value 
agrees with the simulation result described below, and is 
our only fit parameter. 

III. SIMULATION 

Several methods have been proposed in the litera- 
ture to simulate the dense random packing of hard 
spheres. One method used frequently was introduced 
by Lubachevsky and Stillinger,[16[ who simulated hard 
spheres by Molecular Dynamics and slowly compress the 
system until it jams. There is a drawback to this method, 
namely that ever smaller time steps need to be taken as 
close packing is approached. Moreover, the physical re- 
laxation time of the system diverges near the close pack- 
ing density, [l7j and consequently long runs are necessary. 
To circumvent this problem O'Hern et aZ.[18[ used soft 
spheres, and located the minimum energy by a conjugate 
gradient (CG) algorithm. This method is much faster 
than a hard sphere simulation, but the disadvantage is 
that the CG algorithm only simulates the high friction 
limit. 

In general the dense random packing dependends on 
friction, 6] and dissipative interactions may play a role 
not just at hard sphere contact. Granular particles sus- 
pended in a viscous medium also dissipate energy via 
long-range hydrodynamic interactions. Hence we antici- 
pate that the dense random packing also depends upon 
solvent viscosity and on the range of friction. Therefore 
we introduce a new simulation method to include these 
effects. 

Following O'Hern et aL[l8[ and Groot and 
Stoyanov[19| we simulate repulsive elastic spheres 
in the limit T — > 0. The generalization of the repulsive 
force in this model to elastic spheres of unequal size is 

F^ = 2ER ij (D ij -r)6(D ij -r). (4) 

Here, r is the distance between particle centres, Dij = 
(Di+Dj)/2 is the mean diameter and Rij is the harmonic 
mean radius given by Rij — ^DiDj j Dij . The parameter 
E is proportional to the linear elastic modulus of the 
particles [l9|. Henceforth we use E — 1000. 

Instead of using a CG algorithm to search the energy 
minimum, or imposing energy dissipation at particle col- 
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FIG. 3: Volume fraction obtained by extrapolating the P — <f> 
curve to zero pressure. 

lision, we introduce a soft friction function of finite range 
that represents the hydrodynamic interaction between 
spheres. A soft friction has been introduced before to 
simulate hydrodynamics in fluids, in the context of Dis- 
sipative Particle Dynamics [20l l2ll |22| . However, appli- 
cation to particles of unequal size is new to our knowl- 
edge, and because the dense random packing depends on 
friction some care must be taken in defining the friction 
function. 

The most general distance-dependent friction is 

VIP* = -7«v£0(r/r c ) (5) 

where 7^ is a friction factor that may depend on both 
particle sizes, is the radial velocity difference, g(r/r c ) 
is a distance dependent function, and r c is a cut-off dis- 
tance that may again depend on particle size. 

To demonstrate the importance of the friction func- 
tion, we first concentrate on monodisperse systems and 
study the influence of friction range and strength, and 
of particle size. We use a periodic 10 x 10 x 10 box, 
containing from 1222 up to 1290 particles. All par- 
ticles have diameter 1 and mass 1, and interact with 
a repulsive force F — 10 3 (1 — r) for r < 1. First 
we study the influence of the range of the friction in- 
teraction. To this end we use the friction interaction 
F frid = _ v ^(i _ r/r c ) 2 /(l - l/r c ) 2 ; and the range is 
varied as r c = 1.5, 2, 3. With this choice for the fric- 
tion force, the friction at particle contact is unity for all 
systems. All systems are evolved over 10 4 steps or more, 
with St = 0.01. For each system the pressure at T = 
was averaged over 5 independent starting configurations. 

Even though all systems have the same friction 
strength at r — 1, the mere range of friction appears 
to influence the pressure in the final state. As the fric- 
tion range increases, so does the pressure at T = 0. In 
particular, the volume fraction to which the pressure ex- 
trapolates to zero - the dense random packing - varies 



systematically with the force range, see Fig. [3] Although 
the effect is not very large (about 1% variation), it is clear 
that the friction range does have an influence. The ex- 
trapolated value to r c = D = 1, <j> cp = 0.6392 ± 0.0004, 
compares well with the reported mean-field result for 
hard spheres with a friction interaction at contact, Q 
4*rcp — 0.634. Thus we conclude that, to eliminate the 
influence of the friction range, the range must be scaled 
relative to the particle size. 

Next, to demonstrate the influence of particle size and 
friction strength, a series of simulations was done where 
the ratio of the friction range relative to the particle di- 
ameter was kept fixed at r c /D = 1.5. The friction force 
used in this study was F^" c * = — 7/Vy(l — r/r c ) 2 , where 
7/ is a fixed friction factor, to be specified as input vari- 
able. Two sizes were studied, D = 1/2 and D = 1, and 
two friction factors were used, 7/ = 1 and 7/ = 4. In 
these simulations all particle masses were put at m = 1. 
The box size was taken as V = 8 3 , and the conservative 
force was taken as F — 10 3 D(D — r) for r < D, the same 
as in the previous simulations. 

The results are shown in Fig. The red lines give the 
results of low viscosity (7/ = 1) and the black lines give 
the results of higher viscosity (7/ = 4) . Results for small 
particles are denoted by open symbols and dashed curves, 
while results for big particles are denoted by closed sym- 
bols and full curves. This shows that both the friction 
factor and particle size influence the pressure in the glassy 
state. Consequently, the dense random packing density 
must depend on particle size. Even though the effect is 
small it is important, as it points at a reason why the 
dense random packing density is ill-defined. For prac- 
tical reasons we wish to define a dense random packing 
density that does not depend on particle size, i.e. that 
is scale invariant. To obtain this, it is not sufficient to 
have a friction range that is proportional to the particle 
diameter; the packing density depends in a complicated 
way on particle size and on the strength and range of the 
friction force. Empirically there may seem to be some 
scaling when the friction is increased with the square of 
particle size (results for D = 1/2, 7/ = 1 in Fig. 0^ par- 
tially overlap with D = 1, 7/ = 4), but the slopes of the 
curves are clearly different. 

The above results show that we cannot just take any 
friction function. It must reflect the physical properties 
of the hydrodynamic interaction. One physical property 
is the scaling of the hydrodynamic force with particle 
size. On dimensional grounds the friction factor must be 
proportional to particle radius, as in Stokes' law, F = 
6TTi]av. More generally, the (radial) squeeze mode of the 
hydrodynamic force between two rigid spheres at close 
contact behaves asj23jj 

VIP* = -7V^#/fly) « '^<rf ( 6 ) 

where 7 is the friction factor that is proportional to fluid 
viscosity, and h = r — Dm is the distance of closest ap- 
proach. The function g(x) ~ (l/ x ) i s a scaling function 
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FIG. 4: Mean pressure at T = as function of particle vol- 
ume fraction; a (top), for big particles (full symbols) and small 
particles (open symbols) and for low (red) and high (black) 
viscosity; b (bottom) big particles (red dots) and small par- 
ticles (black circles) using friction as in Eq. (0 and mass 
m = D. 



that represents the lubrication force. 

There is a large body of evidence showing that cor- 
rect long-range inertial hydrodynamics is generated even 
if the (divergent) lubrication force between particles is 
replaced by a finite distance-dependent friction[2l|, [22l |. 
It is important however to choose the harmonic mean ra- 
dius as the scaling length (unlike the conjecture by Kim 
and Karilla 23]), otherwise the friction between very un- 
equal spheres would vanish if we remove the divergence 
of g{x). Thus, we use the friction function 



? frict _ 



-7Vy-Rii(l 



h/R t] ) 2 8(R tJ - h). (7) 



which captures the right physics regarding the scaling of 
the range and strength of the viscous interaction with 
particle size. 

Now we can turn to the problem of defining a size 
invariant dense random packing density. Therefore we 
analyse the relevant time scales of the problem. For a 
monodisperse system of elastic particles, a first time scale 



is the oscillation time, t e i = 27r(m/ ED) 1 / 2 . This is the 
elastic time scale. The second time scale in the system is 
the drag relaxation time, td — m/{^D). The dense ran- 
dom packing can only be independent of particle size if we 
maintain a constant ratio between these two time scales, 
t e i/t d = 2nj(D/Em) 1 / 2 . Thus, we have scale invariance 
only if the friction factor satisfies 7 oc (m/D) 1 / 2 . Since 
for most systems mass scales as m oc D 3 , this implies that 
for (soft) elastic spheres with hydrodynamic interaction 
the dense random packing fraction is (weakly) particle 
size dependent. 

To obtain a well-defined dense random packing we are 
forced to choose the particle mass m proportional to D. 
When this choice is made, the above ratio of time scales 
becomes particle size independent, and consequently the 
dense random packing fraction is well-defined. This 
choice has been made henceforth. The predicted scal- 
ing was checked by simulation and holds exactly. The 
pressure as function of time for systems of particle di- 
ameter D = 1/2 and D = 1 fall on top of each other if 
we scale the friction range with particle size (as in Eq. 
([7])) and simultaneously impose m oc D. To demonstrate 
the improvement in system definition, two series of sim- 
ulations were done, again for monodisperse systems of 
particle diameter D — 1/2 and D = 1, with repulsive 
force F — 10 3 D(D — r) for r < D. For a realistic hy- 
drodynamic scaling we used Eq. (J7J), with friction factor 
7 = 0.74. Because the systems are monodisperse the 
cut-off distance for the friction interaction is r c = 1.5-D, 
as in the previous case. To obtain scale invariance, we 
choose the masses as m = D. The systems had fixed 
volume V = 10 3 for D = 1 and V = 5 3 for D = 1/2. 
The systems were integrated over 5000 steps with time 
step St = 0.01 and the pressure was averaged over 10 in- 
dependent runs. The comparison between small and big 
particles is shown in Fig. \Qp. The predicted scaling is 
followed excellently. 

Now that the system is well-defined, we can define a 
fast algorithm to obtain the close packing density. We use 
a variation of the Lubachevsky-Stillinger algorithm [l6j. 
where we make use of the relative softness of the inter- 
action potential. We prepare the system in a random 
conformation (with particle overlaps) and then evolve it 
in an (N, V, T) ensemble until we have a completely equi- 
librated state. For the parameters 7=1 and St = 0.01 
that we used, this requires 3 — 50 x 10 3 time steps. Then 
we switch to an (N, P, T) ensemble, where the pressure is 
steered towards P — 0.01, which is close enough in prac- 
tice to P = (the error in max is of the order 10 -5 ). 
If during a run the pressure falls below 0.001 we switch 
to the L-S algorithm and compress the system in small 
steps until the pressure turns positive. The advantage of 
this method over the standard L-S algorithm is that the 
(high) pressure in the initial (N, V, T) simulation quickly 
drives the system towards P = 0. The final (N, P, T) 
simulation serves to run down the P — <p curve (see Fig. 
QJ)) to locate the intercept at P — 0. Some minor evo- 
lution can however still be observed at P = 0.01. To 
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FIG. 5: a (top): Close packing density as function of the 
friction factor 7. Each data point is an average over 10 inde- 
pendent runs; b (bottom) same data, plotted to 



gain further simulation speed we combined a linked cell 
neighbour search with a Verlet neighbour list [13]. In the 
late stages of evolution, when particles hardly move, this 
leads to a large increase in simulation speed, particularly 
for systems of large particle size difference. 

Even though the close packing density is now well- 
defined, it still depends on friction, or the cooling rate, 
see Fig. [5] In fact this is the very source of the previously 
found dependence of the close packing density on parti- 
cle size and mass. To study this relation, monodisperse 
systems of 6000 particles were used. All particles have di- 
ameter 1 and mass 1, and interact with a repulsive force 
with E = 10 3 . We insert the particles in a box of size 
16.83 (or <f> = 0.66), pre-equilibrate for 5 x 10 4 to 2 x 10 5 
time steps of 5t — 0.01 until the pressure has fully equili- 
brated. Then we run a constant pressure ensemble, steer- 
ing the pressure to P = 0.01, until the volume fraction 
is stable over four decimal places (1.5 x 10 5 time steps). 
All results are averaged over 10 runs. Polycrystallinc 
domains only start to occur for 7 < 0.03; all systems 
shown here are isotropic. Over about a decade we find a 
log(7) dependence for 7 — > (see Fig. [Sk.), and over two 



decades we find a power law decay <j> ~ 0.64 + 0.0028/^/7 
for 7 > 1/4 (see Fig. [5Jd). Therefore we have to make a 
choice for the friction factor, and only refer to the pack- 
ing fraction at that value of 7. Our default value used in 
the next section is 7 = 1. 

For comparison we also evolved a system from its ini- 
tial conformation to equilibrium using a steepest descent 
method in an [N, V, T) ensemble, which should com- 
pare well with the CG algorithm [lg. The final pres- 
sure coincides with our result at 7 — > 00, which demon- 
strates that CG simulates the high viscosity limit. For 
truly hard spheres the modulus diverges, hence the ra- 
tio tei/td = 2tt7(D / 'Em) 1 / 2 -> 0. To simulate this with 
particles of finite modulus, low values of 7 would be pre- 
ferred. This implies that the present method is closer to 
the physical case than the CG algorithm to generate hard 
spheres conformations, unless largely inelastic collisions 
are pertinent. 



IV. RESULTS 

The simple approximation for T described in section 
II is now compared with the results from the sphere 
packing simulation method of section III. Consider first 
bidisperse spheres, as studied for example by Clarke and 
Wiley [25l| and by Yerazunis et al.[2^|, where the larger 
spheres have R times the radius of the smaller, so that 
P 3D (D) cx R 3 (l - w)8{D - l/R) + wS{D - 1). The sim- 
ulation results for binary mixtures are shown in Table [H 
and in Fig. [5J together with the theoretical prediction. 
The big particles have diameter Di = 1, and the small 
particle diameters are D2 = 0.5, 0.3, 0.2 and 0.1. The 
dash-dot curve gives the exact upper limit of the volume 
fraction, Eq. j3|. For diameter D2 > 0.2 we used 6000 
particles; for D2 = 0.1 we used up to TV = 49950 particles 
(at w — 0.8) to prevent finite size effects. All runs were 
evolved over a minimum of 150 000 time steps, and con- 
vergence of the volume fraction was checked by extending 
the evolution of selected systems to 450 000 steps. All 
volume fraction results shown in Fig. [5] are stable up to 
four decimal places. 

A bidisperse system, with moderate to large size ratios, 
has two distinct regimes (and a non-trivial crossover be- 
tween them): When the proportion of large spheres is 
low, they are isolated from one another, like holes in a 
Swiss cheese; while the small spheres form a close-packed 
phase (the 'cheese') between them. On the other hand, 
when the proportion of large spheres is high, these form a 
close-packed structure, leaving the small spheres to 'rat- 
tle around' in the gaps between them. Recent theories 
of the close packed state of bidisperse spheres [TJ, [27j do 
not address the 'rattler' regime adequately. In contrast, 
our theory captures both regimes (exactly, in the limit of 
infinite size ratio), and also the analogous regimes which 
are produced for larger numbers of size classes, such as 
tridiperse spheres. 

Next, we consider a tridisperse distribution. In partic- 
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FIG. 6: Maximum packing fraction for bidisperse spheres of 
different size. R is the size ratio, and w = <^iar g c/(^>iarg C + 
<?Waii) is the relative volume fraction of the large spheres. 
Symbols are simulation results, and solid lines are theoretical 
predictions, based on 20000 rods. The dashed curves give the 
upper limit for infinite size ratio, Eq. (|3J). 



TABLE I: Simulation results for the maximum packing frac- 
tion of bidisperse sphers, with diameters D\ and D^- The 
mass fraction present in the large spheres is given by w — 
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FIG. 7: Maximum packing fraction for a tridisperse distri- 
bution of spheres with size ratios 1:3:9. The composition di- 
agram (a, top) is based on weight fractions; contour lines 
connect points of equal volume fraction, with bold lines at 
volume fractions of 0.65, 0.7, 0.75 and 0.8. The compositions 
used in the simulations are marked by the circles; b (bot- 
tom) shows the comparison between theory and simulation 
per sample point. 



ular we consider the case where the sphere diameters are 
in the ratio 1:3:9. Fig. [7] shows the theoretical predic- 
tion compared with simulation results for selected points. 
The number of particles chosen varies from 6000 to 13400. 
Particular care was taken at large weight fractions of big 
particles, where these may form a stress bearing network. 
A finite size study showed that in that case (specifically 
sample point 1) the number of large particles in the sys- 
tem needs to be above 175 for reliable results. For sample 
point 1 we used 209 big particles from a total of 13300. 
Again, the theory compares very well with the simula- 
tion results. The differences may well be attributed to 
remaining (minor) finite size effects. 

Finally, we consider a log-normal distribution, which is 
defined as P 3D (D) cx exp{-[\n(D/D )] 2 /2a 2 )/D. Thus, 
from Eq. ([2]), we find the rod distribution in the theory 
as P\d{L) oc L erfc[ln(L/£>o)/o'V / 2]. In these simula- 



tions we set an upper diameter to the particles D = 1, 
and choose the mean diameter such that less than 0.1% 
of the particles exceed this size. We used 6000 particles 
that were evolved over 50 000 to 200 000 steps, until the 
volume fraction had converged up to 4 decimal places. 
Fig. [5] shows the comparison between theory and simu- 
lation data. The system at cr = 0.6 is shown in Fig. [T] 
Again, the theory compares very well with the simulation 
results. 



V. CONCLUSIONS 

Concluding, we have introduced a theory for the close 
packing density of hard spheres of arbitrary size distribu- 
tion, based on a mapping onto a one-dimensional prob- 
lem. To test the theory we simulated the dense random 
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FIG. 8: Maximum packing fraction for a log-normal distribu- 
tion. The spread in log(radius) is denoted by a. Simulation 
results are given by symbols; the solid line is the theoretical 
prediction, based on 20000 rods. 

packing of (soft) elastic spheres with hydrodynamic fric- 
tion in 3D in the limit T — > 0, which approaches the hard 
sphere system. For the distributions studied we obtain 
excellent agreement between theory and simulation. The 
theory reproduces the infinite size ratio limit for bidis- 
perse spheres. Hence we expect that this approxima- 
tion will prove useful for more general size distributions. 



The simple structure of the approximation may also be 
amenable to further analysis, and open up new avenues 
for analytical approximations. 

Application of this theory to other space dimensions 
than 3D is straightforward. However, a comparison be- 
tween theory and simulations showed that, although the 
theory is qualitatively correct in 2D, it does not repro- 
duce the simulations as accurately as in 3D. This may 
be related to the mean-field character of the theory, i.e. 
the explicit spatial correlation is lost in the theory. We 
therefore speculate that the theory will be accurate in 3D 
and in higher dimensions. 

The simulations show a weak dependence of the dense 
random packing on fluid viscosity, if the friction force is 
of hydrodynamic origin. In general the dense random 
packing density also depends on particle size, mass and 
elastic modulus. For particles of diameter D, mass m 
and elastic modulus E, suspended in a liquid of viscosity 
rj, we infer that the dense random packing density should 
be a function of the dimensionless group Q = rj 2 D/ (Em). 
For systems of the same Q value but different size, mass 
and viscosity we find excellent scaling of the pressure as 
function of time. Therefore we conclude that the size 
dependence of dense random packing is a kinetic effect 
that disappears when m cx D. Although such scaling is 
artificial, it leads to a well-defined dense random packing, 
which is prudent to test theories that are based only on 
geometrical considerations. 
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